Shocks in nonlocal media 
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We investigate the formation of collisionless shocks along the spatial profile of a gaussian laser 
beam propagating in nonlocal nonlinear media. For defocusing nonlinearity the shock survives 
the smoothing effect of the nonlocal response, though its dynamics is qualitatively affected by the 
latter, whereas for focusing nonlinearity it dominates over filament at ion. The patterns observed in 
a thermal defocusing medium are interpreted in the framework of our theory. 



Shock waves are a general phenomenon thoroughly in- 
vestigated in disparate area of physics (fluids and water 
waves, plasma physics, gas dynamics, sound propagation, 
physics of explosions, etc.), entailing the propagation of 
discontinuous solutions typical of hyperbolic PDE mod- 
els They are also expected in (non- hyperbolic) 
universal models for dispersive nonlinear media, such as 
the Korteweg-De Vries (KdV) and nonlinear Schrodinger 
(NLS, or analogous Gross-Pitaevskii) equations, since hy- 
drodynamical approximations of such models hold true 
in certain regimes (typically, in the weakly dispersive or 
strongly nonlinear case) 0, 0, 0- However, in the lat- 
ter cases, no true discontinuous solutions are permitted. 
The general scenario, first investigated by Gurevich and 
Pitaevskii [3], is that dispersion regularizes the shock, de- 
termining the onset of oscillations that appear near wave- 
breaking points and expand afterwards. This so-called 
collisionless shock has been observed for example in ion- 
acoustic waves [6], or wave-breaking of optical pulses in 
a normally dispersive fiber [7], and recently in a Bose- 
Einstein condensate with positive scattering length 

In this Letter we investigate how nonlocality of the 
nonlinear response affects the formation of a collisionless 
shock in a system ruled by a NLS model. In fact nonlocal- 
ity plays a key role in many physical systems due to trans- 
port phenomena and finite range interactions (e.g. as in 
Bose-Einstein condensation), and can be naively thought 
to smooth and eventually wipe out steep fronts character- 
istic of shocks. More specifically, we place this problem 
in the context of nonlinear optics where nonlocality arises 
quite naturally in different media @, Ell EH, EH , study- 
ing the spatial propagation of a fundamental (gaussian 
TEMoo) laser mode subject to diffraction and nonlocal 
focusing/defocusing action (Kerr effect). In a defocus- 
ing and ideal (local and lossless) medium, high intensity 
portions of the beam diffract more rapidly than the tails 
leading, at sufficiently high powers, to overtaking and os- 
cillatory wave-breaking similar (in ID) to what observed 
in the temporal case [18]. We find that, while shock is 
not hampered by the presence of (even strong) nonlocal- 
ity, the mechanism of its formation as well as post-shock 



patterns are qualitatively affected by the nonlocality. Ex- 
perimental results obtained with a thermal defocusing 
nonlinearity are consistent with our theory and shed new 
light on the interpretation of the thermal lensing phe- 
nomenon. 

Importantly, our theory permits also to establish that 
nonlocality allows the shock to form also in the focusing 
regime where, contrary to the local case, it can prevails 
over filamentation or modulational instability (MI). 

Theory We start from the paraxial wave equation 
obeyed by the envelope A of a monochromatic field 
E = (-^)^ 2 Aexp(ikZ - iuT) (\A\ 2 is the intensity) 
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where k = k$n = is the wave-number, and 
the intensity loss rate. A sufficiently general nonlocal 
model can be obtained by coupling Eq. (1) to an equa- 
tion that rules the refractive index change An of nonlin- 
ear origin. Introducing the scaled coordinates x,y,z = 
X/wo, Y/wo, Z/L, and complex variables i/j = A/y/lo 
and 6 = koL n iAn, where L n \ = (&o|^2|^o) _1 is the non- 
linear length scale associated with peak intensity Iq and 
a local Kerr coefficient ri2 (An = n2\A\ 2 ), Ld = kw^ is 
the characteristic diffraction length associated with the 
input spot-size wo, and L = \JL n \Ld, such model can be 
conveniently written as follows [12[ 



dz 



7 \_^ + xH = - 

-o- 2 V\0 + 6 



.a 



(2) 
(3) 



where a = a L, = d 2 + d 2 , x = ^2/^2 1 = =•=! is 
the sign of the nonlinearity, and a 2 is a free parameter 
that measures the degree of nonlocality. The peculiar 
dimensionless form of Eqs. ([2j[3]) where e = L n \/L = 
yjL n il 'Ld is a small quantity, highlights the fact that we 
will deal with the weakly diffracting (or strongly non- 
linear) regime, such that the local a = and lossless 
a = limit yields a semiclassical Schrodinger equation 
with cubic potential (s and z replace Planck constant 
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and time, respectively). We study Eqs. ([2]|3]) subject 
to the axi- symmetric gaussian input ipo(r) = exp(— r 2 ), 
r = \/x 2 + y 2 , describing a fundamental laser mode 
at its waist. For e <C 1, its evolution can be studied 
in the framework of the WKB trasformation ^(r, z) = 
y/p(r, z) exp [i<j>(r, z) j e] Q. Substituting in Eqs. d2H3j) 
and retaining only leading orders in £, we obtain 
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= -ap; u z + uu r - x&r = 0, 
9 = p. (4) 



where u = (f) r is the phase chirp, and D = 2 is the trans- 
verse dimensionality. The ID case described by Eqs. (|4]) 
with D = 1 and r — > x (<9 y = 0) illustrates the ba- 
sic physics with least complexity. In the defocusing case 
(x = —1) for an ideal medium (a = a = 0, 9 = p), 
Eqs. (|4j) are a well known hyperbolic system of conser- 
vation laws (Eulero and continuity equations) with real 
celerities (or eigenspeeds, i.e. velocities of Riemann in- 
variants) v ± = u ± \J~XPi which rules gas dynamics [u 
and p are velocity and mass density of a gas with pres- 
sure oc p 2 ). A gaussian input is known to develop two 
symmetric shocks at finite z 0] . Importantly the diffrac- 
tion, which is initially of order £ 2 , starts to play a major 
role in the proximity of the overtaking point, and regu- 
larize the wave-breaking through the appearance of fast 
(wavelength ~ e) oscillations which connect the high and 
low sides of the front and expand outwards (far from the 
beam center) Such oscillations, characteristic of a col- 
lisionless shock, appear simultaneously in intensity and 
phase chirp (u) as clearly shown in Fig. QJa,c). 

In the nonlocal case, the index change 9{x) can be 
wider than the gaussian mode (for large a) and the shock 
dynamics is essentially driven by the chirp u with p adia- 
batically following. This can be seen by means of the fol- 
lowing approximate solution of Eqs. (|4]): considering that 
the equation for p is of lesser order [0(e)], with respect 
to those for 9 and u [0(1)], we assume p = exp(— 2x 2 ) 
unchanged in z and solve exactly the third of Eq. (|4|) for 
9(x) (though derived easily, its full expression is quite 
cumbersome). Then, applying the theory of characteris- 
tics the second of Eqs. ([4]) is reduced to the following 
ODEs, where dot stands for d/dz 



u; u = x®x, 



(5) 



equivalent to the motion of a unit mass in the potential 

( 

V(x) = ~x9 with conserved energy 
The solution of Eqs. ([5|) with initial condition x(0) = 
s, u(0) = yields x(s, z), u(s, z) in parametric form, from 
which overtaking is found whenever u(x, z) (obtained 
by eliminating s) becomes a multivalued function of x 
at finite z = z s . The shock point corresponding to 
\du/dx\ — > oo is found from the solution u(x, z) displayed 
in Fig. [2ja) [121(b)], at positions x = ±x s ^ (defocusing 
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FIG. 1: (Color online) ID spatial profiles of phase chirp u(x) 
(a-b) and amplitude |^(x)| = sj p(x) (c-d), as obtained from 
Eqs. (12113)) with e = 10" 3 , a = 0, x = -1 (defocusing), and 
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for different z as indicated: (a-c) local case, 



a = 0; (b-d) nonlocal case, a — 5. 




FIG. 2: (Color online) (a) u(x) for different z and x — 1 
(focusing), a = 1; (b) as in (a) for x — — 1 (defocusing); (c) 
shock distance z s (x = — 1 bold solid, x — 1 thin solid) and 
shock position x s in the defocusing case (dashed line). 



case) or x s = (focusing case). The shock distance z s 
increases with a in both cases, as shown in Fig. [21(c). 

We have tested these predictions by integrating nu- 
merically Eqs. d2H3|). Simulations with x — — 1 [ see Fig. 
[U(b,d)] show indeed steepening and post-shock oscilla- 
tions in the spatial chirp u, which are accompanied by 
a steep front in p moving outward. The shock location 
in x and z is in good agreement with the results of our 
approximate analysis summarized in Fig. [21 

Numerical simulations of Eqs. d2E|) validates also 
the focusing scenario. The field evolution displayed in 
Fig. [31(a) exhibits shock formation at the focus point 
(x s = 0, z s ^8, for a = 5) driven the phase whose chirp 
is shown in Fig. (3](b). This is remarkable because, in 
the local limit a = 0, the celerities become imaginary 
(the equivalent gas would have pressure decreasing with 
increasing density p), and no shock could be claimed to 
exist. In this limit, the reduced problem (|4j) is elliptic and 
the initial value problem is ill-posed [13], an ultimate con- 
sequence of the onset of MI: modes with transverse (nor- 
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FIG. 3: (Color online) Level plot of the intensity in the 
focusing case (x — 1, s — 0.01): (a) nonlocal case (a 2 — 
25); (b) chirp profile for various z for (a); (c) quasi- local case 
(a 2 = 10- 5 ). 



malized) wavenumber q < Aq grow exponentially with z, 
with both gain g and bandwidth Aq scaling as 1/e. How- 
ever, the nonlocal response tends to frustrate MI (see also 
Refs. d, HH), as shown by standard linear stability anal- 
ysis which yields g = ^/ d(2x — d)/e 2 (we set d = e 2 q 2 /2 
and x = x/(l + ^Q 2 ))^ m turn implying a strong reduc- 
tion of both gain and bandwidth for large a. In order 
to emphasize the difference between the local and non- 
local regime, we contrast Fig. [3ja) with the analogous 
evolution [see Fig. 09(c)] obtained in the quasi-local limit 
(a 2 = 10 -5 ), which appears to be clearly dominated by 
filamentation. 

Thermal nonlinearity The physics of the defocusing 
case can be experimentally tested by exploiting ther- 
mal nonlinear it ies of strongly absorptive bulk samples, 
that we show below to fit our model. In this case, the 
system relaxes to a steady-state refractive index change 
An = (dn/dT)AT, where dn/dT is the thermal coef- 
ficient, and AT the local temperature change due to 
optical absorption. It is well known that this so-called 
thermal lens distorts a laser beam propagating in the 
medium LLi 
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|16| . However, only perturbative ap- 
proaches to the problem have been proposed (ray optics 
or Fresnel diffraction theory is applied after the lens pro- 
file is worked out from gaussian ansatz [3), while the 
role of shock phenomena was completely overlooked. 
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FIG. 4: 

a = 1: 



2D evolution according to Eqs. (|2]|3) with a = 1, 
(a) radial phase chirp at different z, as indicated, 



showing steepening and shock formation for e — 10 ~ ; (b) cor- 
responding intensity profile \ip\ 2 (maximum scaled to unity) 
at z = 4.9; (c) transverse intensity profile vs. x (at y = 0) at 
z = l/(4e) and different values of e (ao = 62cm -1 , a = 0.3). 



We assume that the temperature field AT = 



ATj_ (X, Y) obeys the following 2D heat equation 

(d 2 x + d\ )AT ± - CAT ± = - 7 \A\ 2 (6) 

where the source term account for absorption pro- 
portional to intensity through the coefficient 7 = 
cto/ipoCpDr), where po the material density, c p the spe- 
cific heat at constant pressure, and Dt is the thermal 
diffusivity (see e.g. [l6[). Eq. (|6]) has been already em- 
ployed to model a refractive index of thermal origin in 
Ref. [10], and in Ref. [Tl| in the limit C = which 
is equivalent to consider the range of nonlocality (mea- 
sured by 1/C, see below) to be infinite. Starting from 
the 3D heat equation V 2 AT = — 7|^4| 2 , the latter regime 
amounts to assume AT(X, Y, Z) = AXj_(X, Y), which is 
justified when longitudinal changes in intensity \A\ 2 are 
negligible as for solitary (invariant in Z) wave-packets 
in the presence of low absorption [11]. Viceversa, in the 
regime of strong absorption, we need to account for lon- 
gitudinal temperature profiles that are known from solu- 
tions of the 3D heat equations to be peaked at charac- 
teristic distance Z in the middle of sample and decay to 
room temperature on the facets [141 ] . Since highly non- 



linear phenomena occurs in the neighborhood of Z where 
the index change is maximum, we can use a (longitudi- 
nal) parabolic approximation with characteristic width 
L e ff(~ L) of the 3D temperature field AT(X, Y,Z) = 

A2\(X, Y) and consequently approximate 
)AT± - L~f f AT±, so that the 3D 
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heat equation reduces to Eq. (|6|) with C = 1/L 2 ^j. 
Following this approach, Eq. (j6]) coupled to Eq. (Q]) 
can be casted in the form of Eqs. ([2]|3]) by posing 
= k L n i\dn/dT\AT ± and a 2 = l/(Cw§) = L 2 eff /w%. 
The model reproduces the infinite range nonlocality for 
negligible losses (L e ff — > 00); while for thin samples 
[\{d 2 x + d$r)AT±\ « \L~f f AT ± \], L eff can be related 
to the Kerr coefficient n 2 as 



L eff = 



n 2 



j\dn/dT\ 



l D T poc p \n 2 \ 
ao\dn/ dT \ 



(7) 



which establishes a link between the degree of nonlocal- 
ity and the strength of the nonlinear response (similarly 
to other nonlocal materials [Hj]). 

Having retrieved the model Eqs. ([2]|3j), let us show 
next that the scenario illustrated previously applies sub- 
stantially unchanged in bulk (2D case) even on account 
for the optical power loss (a 7^ 0). An example of the 
general dynamics is shown in Fig. 0J where we report a 
simulation of the full model ([2](3j), with a 2 = 1 and rela- 
tively large loss a = 1. In analogy to the ID case, Fig. 
Ufa) clearly shows that the radial chirp u = cV steep- 
ens and then develop characteristic oscillations after the 
shock point (z c± 6, where \d r u\ — > 00). Correspond- 
ingly the intensity exhibits also an external front which 
is connected to a flat central region with a characteris- 
tic overshoot [see Fig. HJ^b)] corresponding to a brighter 



ring [inset in Fig. Etc)]. For larger distances this struc- 
ture moves outward following the motion of the shock. 
In the experiment such motion can be observed, at fixed 
physical lenght, by increasing the power, which amounts 
to decrease e while scaling z and a accordingly (z oc 
aoce), as displayed in Fig. [H^c) for a = 0.3. 

As a sample of a strongly absorbing medium we choose 
a 1 mm long cell filled with an acqueous solution of 
Rhodamine B (0.6 mM concentration). Our measure- 
ments of the linear and nonlinear properties of the sam- 
ple performed by means of the Z-scan technique gives 
data consistent with the literature jTjj], and allows us 
to extrapolate at the operating vacuum wavelength of 
532 nm, a linear index n = 1.3, a defocusing nonlin- 
ear index n<i = 7 x 10 -7 cm 2 W _1 , and ao = 62 cm -1 . 
For our sample Dt = 1.5 x 10 -7 m 2 s _1 , po = 10 3 kg 
m- 3 , c p = 4 x 10 s Jkg^K' 1 and \dn/dT\ = 10~ 4 K" 1 
(7 ^ 10 4 K W" 1 ), and exploiting Eq. © we estimate 
L e ff = 10/xm (L e ff « L because of the strong ab- 
sorption that causes strong heating of our sample near 
the input facet), and correspondingly the degree of non- 
locality a = 0.3. We operate with an input gaussian 
beam with fixed intensity waist wqj = wq/ y/2 = 20 /xm 
(Ld = 12 mm) focused onto the input face of the cell. 
With these numbers, an input power P = ttWqjIq = 200 
mW yields a nonlinear length L n \ = 8 fim (L = 0.3 mm), 
which allows us to work in the semiclassical regime with 
e = 0.025. The radial intensity profiles together with the 
2D patterns imaged by means of a 40 x microscope objec- 
tive and a recording CCD camera are reported in Fig. [5l 
As shown the beam exhibits the formation of the bright 
ring whose external front moves outward with increasing 
power, consistently with the reported simulations. We 
point out that, at higher powers, we observe (both ex- 
perimentally and numerically) that the moving intensity 
front leaves behind damped oscillations that correspond 
to inner rings of lesser brightness, as reported in litera- 
ture [15]. This, however, occurs well beyond the shock 
point that we have characterized so far. 

In summary, the evolution of a gaussian beam in the 
strong nonlinear regime is characterized by occurence of 
collisionless (i.e., regularized by diffraction) shocks that 
survive the smoothing effect of (even strong) nonlocality. 
While experimental results support the theoretical sce- 
nario in the defocusing case, the remarkable result that 
the nonlocality favours shock dynamics over filamenta- 
tion requires future investigation. 
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responding 2D output patterns. 
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